# TODO: Add comment
# 
# Author: wuping
###############################################################################
library(FNN)
source("./sourceAll.R")
Nreal=50000
population=new("population",x=runif(Nreal,0,1000),y=runif(Nreal,0,500),plotdim=c(1000,500))
po=ppp(x=population@x,y=population@y, window=owin(c(0,1000),c(0,500)))
po=unique(po)
population@x=po$x
population@y=po$y

result=data.frame("cp"=NA,"q"=NA,"k"=NA,"Nmean"=NA,"Nsd"=NA,"mse"=NA)
for (cp in c(0,1)){
	for (q in c(1,2,4)){
		for (k in c(1,2,3)){
			angleorder=new("angleOrder",n=200/q,q=q,k=k,edge=80,stopLimit=1,composite=as.logical(cp))
			N=numeric()
			for(i in 1:2000){
				N[i]=pSizeTor(sampling(angleorder,population,quick=T))
			}
			N=N[!is.na(N)]
			mse=sqrt(sum((Nreal-N)^2)/length(N))/Nreal
			tp=c(cp,q,k,mean(N),sd(N),mse)
			
			result=rbind(result,tp)
		}
	}
}

result=result[-1,]
print(result)
print(result[order(result$mse),])

#write.table(result,"PDEs for random distribution.dat",quote=F)



####
# 表现好坏（mse的大小）首先取决于q与k的乘积
# 其次还受边缘效应的影响，并且边缘效应随着k和q的增加其影响也逐步增强。

#当lambda很大，边缘缓冲区远大于期望邻体距离时，应该是k=3，q=4最好。
#同时由于组合指数的方差大于相同k、q下的单一方差，因此在随机分布时
#组合的表现应该会比单纯的要差，
#当这种差可能会在实际的聚集格局分布中补偿回来，使得组合的更好！

#但当固定n*p的乘积时，只有增加k才能最终改善估计结果，而增加q增加了边缘效应的干扰可能性，
#使得估计结果更差

q=c(1:100)
k=c(1:100)
k1=gamma(k+1/2)/gamma(k)
q1=sqrt(q)
data.frame(x1=k1[-1]-k1[-length(k)],
x2=q1[-1]-q1[-length(q)])

#方差比较
v1=1/(nq*k-1)
v2=gamma(nq*k/2)^4/gamma(nq*k/2-1/2)^4/(nq*k/2-1)^2-1
#exp(4*lgamma(nq*k)-4*lgamma(nq*k-1/2)-2*log(nq*k-1))-1
